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finite  rise  time  the  amplitude  of  the  wave  passing  through  the  surface  atoms 
is  diminished  when  compared  with  the  zero-time  case.  For  the  anharmonic 
lattice  the  head  of  the  wave  develops  into  a solitary  wave  train  with  an 
oscillatory  tail,  and  for  certain  rise  times  and  anharmonicity  parameters  an 
envelope  soliton  forms  behind  the  shock  front.  This  envelope  soliton  travels 
much  slower  than  the  shock  wave.  The  relevance  to  Army-related  problems  is 
discussed. 
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1 .  INTRODUCTION 


Historically,  the  most  common  assumption  for  generating  a shock  wave 
in  a lattice  from  a microscopic  point  of  view  is  to  maintain  the  end 
particle  at  a steady  compression  velocity  u for  all  times.  Paskin  and 


1-2 


..  -3-6 

Manvi 


Dienes 
and  Powell  and  Batteh 


Duval 1 
13 


3-5 


and  Lowell 


3-4 


_ .7-12 

Tasi  , 


and  Mus grave 
14 

have  used  the  above  assumption  while  Beckett 


1.1-12 


IT  A.  Paskin  and  G.J.  Dienes,  "Molecular  Dynamic  Simulations  of  Shock 
Waves  in  a Three-Dimensional  Solid",  J.  Appl.  Phys.  £3,  1605 
(1972). 

2.  A.  Paskin  and  G.J.  Dienes,  "A  Model  for  Shock  Waves  in  Solids  and 
Evidence  for  a Thermal  Catastrophe",  Solid  State  Comm.  17_,  197 
(1975). 

3.  R.  Manvi,  G.E.  Duvall,  and  S.C.  Lowell,  "Finite  Amplitude  Longi- 
tudinal Waves  in  Lattices",  Int.  J.  Mech.  Sci.  ££,  1 (1969). 

4.  G.E.  Duvall,  R.  Manvi,  and  S.C.  Lowell,  "Steady  Shock  Profile 
in  a One- Dimensional  Lattice",  J.  Appl.  Phys.  40,  3771  (1969). 

5.  R.  Manvi,  and  G.E.  Duvall,  "Shock  Waves  in  a One- Dimensional, 
Non-Dissipating  Lattice",  Brit.  J.  Appl.  Phys.  2_,  1389  (1969). 

6.  R.  Manvi,  "Shock  Wave  Propagation  in  a Dissipating  Lattice  Model", 
Ph.D.  Thesis  (Washington  State  University,  1968)  (Unpublished). 

7.  J.  Tasi,  "Perturbation  Solution  for  Growth  of  Nonlinear  Shock 
Waves  in  a Lattice",  J.  Appl.  Phys.  4£,  4016  (1972).  See  also 
Erratum  (J.  Appl.  Phys.  44,  1414  (1973)). 

8.  J.  Tasi,  "Perturbation  Solution  for  Shock  Waves  in  a Dissipative 
Lattice",  J.  Appl.  Phys.  44,  2245  (1973). 

9.  J.  Tasi,  "Far-Field  Analysis  of  Nonlinear  Shock  Waves  in  a Lattice", 
J.  Appl.  Phys.  44.  4569  (1973). 

10.  J.  Tasi,  "Reflection  of  Nonlinear  Shock  Waves  in  a Lattice",  J. 

Appl.  Phys.  47,  5336  (1976). 

11.  M.J.P.  Musgrave  and  J.  Tasi,  "Shock  Waves  in  Diatomic  Chains  - I. 
Linear  Analysis",  J.  Mech.  Phys.  Solids  24_  19  (1976). 

12.  J.  Tasi  and  M.J.P.  Musgrave,  "Shock  Waves  in  Diatomic  Chains  - II. 
Nonlinear  Analysis",  J.  Mech.  Phys.  Solids  24_,  43  (1976). 

13.  J.  Powell  and  J.  Batteh,  "Shock  Propagation  in  the  One-Dimensional 
Lattice",  BRL  Report  No.  2009,  1977.  (AD  #A044791) 

14.  D.H.  Tsai  and  C.W.  Beckett,  "Shock  Wave  Propagation  in  Cubic 
Lattices",  J.  Geophys.  Res.  7£,  2601  (1966). 
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MacDonald1  , and  Tsai1  ” have  assumed  that  it  is  more  realistic  to 
accelerate  the  end  particle  from  zero  to  its  final  velocity  in  a time 
tmax’  Indeed,  it  seems  impossible  for  the  end  particle  to  achieve  its 

final  velocity  instantaneously.  There  should  be  some  effect  on  wave 
propagation  in  the  surface  atoms,  and  perhaps  some  other  phenomena  owing 
to  the  differences  in  the  two  cases.  However,  it  is  not  clear  whether 
the  results  will  be  changed  enough  to  warrant  the  additional  complications 
of  a more  realistic  model.  Therefore,  the  purpose  of  this  report  is  to 
obtain  a better  understanding  of  wave  propagation  in  a lattice  for  the 
two  different  cases.  A fringe  benefit  of  this  study  is  that  some  results 
have  been  obtained  for  the  most  common  assumption,  which  will  be  reported 
here  for  the  first  time. 

The  specific  problem  chosen  for  study  in  this  report  is  a small  part 
of  an  ever  increasing  research  effort  to  describe  shock  propagation  in 
solids  from  a microscopic  point  of  view.  Usually  this  nonlinear  problem 
of  solving  Newton's  second  law  for  the  motion  of  individual  particles 
is  done  numerically  on  the  computer.  Since  this  method  starts  from  such 
a fundamental  equation,  many  physical  effects  can  be  treated  without 
the  usual  simplifying  assumptions  of  continuum  mechanics.  For 
instance,  a well-  known  lattice-dynamical  result  predicts  that  if  a one- 
dimensional lattice  with  linear  interatomic  forces  is  subjected  to  steady 
compression,  the  wave  profile  will  spread  as  it  travels  farther  into 
the  lattice.  This  effect  is  explained  by  the  fact  that  the  normal- 
mode frequencies  have  different  group  velocities.  However,  this  disper- 
sion arising  from  the  discrete  nature  of  the  lattice  is  not  included  in 
the  hydrodynamic  approach.  Furthermore,  when  nonlinear  interatomic  forces 
are  taken  into  account,  the  wave  profile  steepens  and  the  normal  modes  of 
the  crystal  become  coupled.  The  consequences  of  this  coupling  can  best 
be  determined  by  a microscopic  approach. 

Many  physical  effects  are  assumed  in  a continuum  approach  to  have 
specific  characteristics  based  on  intuitive  reasoning.  One  such  assump- 
tion is  that  the  shock  compressed  material  has  yielded  completely, so  that 
the  stresses  may  be  assumed  to  be  hydrostatic.  The  usual  justification 
for  tiiis  statement  is  that  the  high  pressure  in  the  compressed  region 

causes  the  shear-yielding  to  be  complete.  However,  Tsai1^  has  suggested 


15.  D.H.  Tsai,  "An  Atomistic  Theory  of  Shock  Compression  of  a Perfect 
Crystalline  Solid",  in  Accurate  Characterization  of  the  High-Pressure 
Environment , edited  by  E.C.  Lloyd,  Natl.  Bur.  Stds.  Spec.  Publ.  No. 
326  (U.S.  GPO,  Washington,  DC,  1971),  p.  105. 

16.  D.H.  Tsai  and  R.A.  MacDonald,  "Second  Sound  in  a Solid  Under  Shock 
Compression",  J.  Phys.  C,  6^,  L171  (1973). 

17.  H.  Prask,  P.  Kemmey,  S.  Trevino,  D.H.  Tsai,  and  S.  Yip,  "Computer 
Simulation  Studies  of  the  Microscopic  Behavior  of  Shocked  Solids", 
Proceedings  of  the  Conference  on  Mechanisms  of  Explosion  and  Blast 
Waves,  editor  J.  Alster  (JTCG/ALNNO,  Naval  Weapons  Station,  Yorktown, 
VA,  1973,  XVI). 
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that  the  yield  strength  of  the  material  is  probably  much  higher  under 
dynamic  conditions.  This  problem  of  time  - dependent  yielding  under 
the  transient  condition  of  shock  compression  has  not  been  exhaustively 
studied  either  experimentally  or  theoretically  in  the  high-pressure 
regime.  Furthermore,  it  is  not  clear  if  the  steady  state  assumption, 
and  thus  the  Hugoniot  relationships  in  their  usual  form,  apply.  Other 
assumptions  for  continuum  treatments  of  shock  propagation  concern  the 

18  19 

equation  of  state  for  the  solid  and  the  nature  of  viscous  effects  ’ 
Actually,  little  is  known  about  the  equation  of  state,  and  the  origins  of 
viscous  effects  are  not  completely  understood.  In  principle,  by  calcu- 
lating the  positions  and  velocities  of  the  particles  as  a function  of 
time  in  the  microscopic  approach,  the  quantities  such  as  pressure, 
density,  and  temperature  can  be  determined  without  an  equation  of  state 
or  assumed  viscous  effects. 


Finally,  nonlinear  (anharmonic)  lattice  dynamics  has  been  greatly 

20 

influenced  by  the  work  of  Fermi,  Pasta  and  Ulam  (FPU)  . At  the  time  of 
their  paper  it  was  generally  accepted  that  small  nonlinearities  would 
lead  to  equipartition  of  energy.  In  a linear  (harmonic)  system  energy 
deposited  into  a normal  mode  can  never  flow  to  another  normal  mode,  but 
it  was  thought  that  a small  nonlinearity  would  cause  energy  to  flow  from 
one  mode  to  another  until  the  time-averaged  energy  of  each  mode  was  the 
same.  FPU  studied  the  vibration  of  particles  connected  by  nonlinear 
springs  by  using  a computer.  They  found  that  the  system  did  not  approach 
thermal  equilibrium,  but  rather  returned  to  its  original  state  after  a 
recurrence  time.  On  the  other  hand,  the  hydrodynamic  theory  assumes 
that  thermal  equilibrium  exists  behind  the  shock  and  allows  for  only  small 
deviations  from  equilibrium  within  the  front.  A recent  report  by  Powell 
13 

and  Batteh  examines  this  question  in  some  detail,  including  a discussion 
of  the  similarities  and  differences  in  the  results  and  interpretations  of 
earlier  authors. 


The  approach  to  thermal  equilibrium  in  the  compressed  region  is 
important  to  the  theory  of  detonation.  For  instance,  the  Zeldovich-von 

21 

Neumann-Doring  (ZND)  theory  of  detonation  treats  the  shock  front  as 
a mathematical  discontinuity.  The  only  function  of  the  shock  wave  is  to 
provide  the  energy  necessary  to  raise  the  temperature,  density,  and 
pressure  of  the  lattice  to  values  higher  than  in  the  undisturbed  lattice. 
The  condensed  region,  where  chemical  reactions  occur,  is  assumed  to  be 
in  thermodynamic  equilibrium.  However,  numerical  and  perhaps  analytical 


18. 

19. 

20. 

21 


W.  Band,  "Studies  in  the  Theory  of  Shock  Propagation  in  Solids", 

J.  Geophys.  Res.  65_,  695  (I960). 

D. R.  Bland,  "On  Shock  Structure  in  a Solid",  J.  Inst.  Math. 
Applications  1_,  56  (1965). 

E.  Fermi,  J.R.  Pasta,  and  S.M.  Ulam,  "Studies  in  Nonlinear  Problems", 
Los  Alamos  Sci.  Lab.  Rep.  LA-1940,  1955;  also  in  Collected  Works 

of  Enrico  Fermi  (Univ  of  Chicago  Press,  Chicago,  1965),  V.  II,  p.  978. 
B. Lewis  and  G.  von  Elbe,  Combustion,  Flames  and  Explosion  of  Gases 
(Academic  Press,  New  York,  1951),  Chap.  XI. 
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solutions  may  show  that  the  shock  profile  is  not  steady,  or  the  time  to 
re-establish  equilibrium  is  of  the  same  order  as  of  the  time  required 
for  a chemical  reaction.  In  that  case  some  assumptions  used  in  detona- 
tion theory  will  have  to  be  changed. 

2.  THEORY 


This  model  treats  a semi-infinite  chain  of  atoms,  of  mass  m,  which 
interact  pairwise  through  a Morse  potential.  At  equilibrium  the  lattice 

spacing  between  neighboring  atoms  is  a.Q,  the  displacement  of  the  j 


atom  from  its  equilibrium  position  is  given  by  the  coordinate  x.,  the 


,th 


3 


distance  to  the  j atom  from  the  origin  of  the  system  located  at  the 
equilibrium  position  of  the  zeroth  atom  is  r.,  and  the  corresponding 


3 


.th 


distance  to  the  equilibrium  position  of  the  j 
Figure  1).  The  coordinates  obey  the  relation. 


atom  is  r 


Oj 


(see 


r . 


J 


r . 


oj 


+ x . . 
3 


(2.1) 


vo  r_a°_*i 

— ► #vwv^vwv^vvv  • • • 

0 1 2 N-l  N 


Figure  1.  Model  for  simulating  shock  propagation  in  a one-dimensional 
discrete  lattice. 


The  aim  of  this  calculation  is  to  obtain  the  solution  for  the  classical 
equations  of  motion  of  the  atoms  when  the  zeroth  atom  is  subjected  to  an 
acceleration  which  changes  its  velocity  from  zero  to  its  final  value  u 
in  a time  tm  . Thereafter,  the  initial  particle  travels  at  a steady 

compression  velocity  u. 


The  Hamiltonian  for  the  chain  can  be  written  as 

oo 

1 2 

H s 2 m I v . + * (r^  r2,  •••,  rR)  , 


(2.2) 


3=0 


■ th 


where  v.  = dx./dt  is  the  velocity  of  the  j particle  and  9 is  the  total 
3 3 

potential  energy  of  the  lattice.  Newton's  second  law  for  the  j particle 
becomes 

d x. 

m -!1  - F ♦ F ext  , (2.3) 

dt  •’  J 
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3$  th 

where  F.  = - -r — is  the  force  exerted  on  the  j particle  by  the  remaining 
J ox . 

^ 6Xt 

atoms  of  the  lattice  and  F.  is  the  corresponding  external  force.  In 

ext  ^ 

this  investigation,  F^  is  zero  for  all  particles  except  j=o,  when  it 
is  varied  in  such  a manner  that  this  particle  goes  from  zero  velocity 


to  the  steady  compression  velocity  u in  a time  t 


max 


It  is  the 


intention  of  this  report  to  investigate  the  difference  between  the  above 

6Xt 

Fq  and  the  more  usual  one  which  constrains  the  zeroth  particle  to 
move  at  constant  velocity  u for  all  times. 

If  we  assume  that  4>  can  be  obtained  by  summing  over  two  body 
potentials,  i.e., 


<‘>(r01’  r02’ 


r.  . 


■)  = l 

i <j 


(2.4) 


where  r^ . = r^-r.,  then  <f>  can  be  expanded  in  a Taylor  series  about  the 

equilibrium  positions  of  the  relative  displacements,  r^  . = rQ^  - r^^ . 

Furthermore,  if  we  assume  that  the  deviations  of  the  r. . from  their 

ij 

equilibrium  value  are  small,  the  expansion  can  be  truncated  a^'ter  second 
order  terms  such  that 


<t=4>  (• 
o 


ij 


T I 

i,  j=o 


^ . XX. 

IJ  1 J 


(2.5) 


where  4>o  is  a constant  which  will  arbitrarily  be  set  equal  ;o  zero  here- 
after and  where 


(2.6) 


The  first-derivative  term  in  the  Taylor  series  expansion  for  the  potential 
vanishes  since  it  is  the  negative  of  the  force  on  the  j atom  in  the 
equilibrium  configuration,  which  is  zero.  The  potential  in  Eq.  (2.5)  is 
called  the  harmonic  potential.  Finally,  if  we  assume  only  nearest-neighbor 

interactions  such  that  only  the  j+lst  and  j-lst  atom  exert  an  appreciable 

force  on  the  j ^ atom,  we  have 


!>..=-  y (6 . . , - 26.  . + 6.  . .) 
1 J i.J-1  ij  i.J  + l 


(2.7) 


15 
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where  y is  the  force  constant  of  the  "spring"  connecting  successive 
particles  and  6 is  the  Kronecker  5.  Equations  (2. 3) -(2. 7)  then  imply 

that  the  equation  of  motion  of  the  j**1  particle  is  given  in  the  harmonic 
approximation  by 


m 


y(x.  ,-2x.+x.  ,)  + F. 
J+1  J J-l'  J 


ext 


(2.8) 


Equation  (2.8)  is  a linear,  second-order,  differential  equation  and  it 
has  an  exact  analytic  solution  for  certain  external  forces,  which  will 
be  given  later. 

However,  while  Equation  (2.8)  can  generally  be  used  for  calculating 
the  equilibrium  properties  of  a lattice  especially  at  low  temperatures, 
it  is  not  a good  starting  point  to  describe  a shock  wave  in  a solid. 

The  first  reason  is  that  at  the  high  temperatures  present  in  shock  waves 
the  relative  displacements  of  the  atoms  from  equilibrium  are  so  large 
that  higher  order  terms  must  be  retained  in  Equation  (2.5).  Second,  in 
the  harmonic  approximation  the  shock  energy  is  initially  distributed 
among  the  normal  modes  in  a nonequilibrium  fashion  and  there  is  no 
coupling  mechanism  allowing  the  crystal  to  thermalize  after  the  shock 
has  passed.  Furthermore,  the  steepening  of  the  wave  profile  in  a 
compressed  lattice  is  caused  by  the  nonlinear  terms.  Therefore,  the 
shock  wave,  which  results  from  the  steepening,  must  include  nonlinear 
terms. 


In  the  present  research  a Morse  potential  was  used,  and  only  nearest- 
neighbor  interactions  were  assumed.  The  Morse  potential  can  be  written 
as 


00 


i=l 


-a(x.-x.  .)  2 

v i i-l'  , 
e -1 


(2.9) 


where  D,  and  a are  constants  which  are  usually  fit  to  the  experimental 
data. 


3.  THE  HARMONIC  LATTICE 

A harmonic  lattice  cannot  support  a shock  wave  in  the  usual  sense 
for  the  reasons  given  in  the  preceding  section.  However,  the  calcula- 
tion was  performed  to  gain  a better  physical  understanding  of  the 
analogous  nonlinear  case.  In  addition,  the  calculations  for  the 
anharmonic  case  must  reduce  to  the  harmonic  case  as  the  nonlinear  terms 
go  to  zero. 

22 

We  shall  begin  by  presenting  the  solution  of  Schroedinger  for  the 


22.  E.  Schroedinger,  "Zur  Dynamik  Elastisch  Gekoppelter  Punktsysteme" , 
Ann.  Phys.  44,  916  (1914). 
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equations  of  motion  for  the  atoms  of  an  infinite,  one-dimensional, 
harmonic  lattice  for  arbitrary  initial  conditions.  When  the  lattice 
is  in  equilibrium  neighboring  atoms  are  uniformly  separated  a distance 
aQ,  where  the  zeroth  atom  is  located  at  the  origin  (see  Figure  2).  All 

the  external  forces  are  zero.  By  observing  the  symmetry  of  the  problem 
we  can  choose  the  initial  velocities  and  positions  of  the  particles  for 
j < 1 so  that  the  zeroth  particle  travels  to  the  right  in  a prescribed 
manner.  The  result  of  the  above  selection  is  to  change  the  boundary- 
value  problem  for  wave  propagation  in  a semi-infinite  lattice  to  an 
initial-value  problem  for  an  infinite  chain.  We  then  use  this  model  to 
investigate  wave  propagation  for  two  sets  of  boundary  conditions  on  the 
zeroth  particle.  In  the  first  case,  the  zeroth  particle  is  set  in 
motion  at  constant  velocity  u,  and  in  the  second  case  it  is  accelerated 
from  zero  to  its  final  velocity  u in  a time  t after  which  it  moves 
at  constant  velocity  u. 

3.1  General  Solution  of  the  Equations  of  Motion  for  the  Semi-Infinite, 
One-Dimensional,  Harmonic  Chain  with  Boundary  Conditions. 


Consider  a chain  consisting  of  N atoms  (see  Figure  2)  connected 
by  harmonic  springs  of  force  constant  y.  Every  atom  has  mass  m and  is 
labeled  by  index  j 


N-l 

2 


< j< 


N-l 

2 


We  will  assume  for  convenience  that  N is  odd.  It  will  be  made  arbi- 
trarily large  in  the  final  results. 

The  differential  equation  of  motion  for  the  jt*1  atom,  assuming 
only  nearest-neighbor  interactions,  is  given  by  Eq.  (2.8)  without  the 
forcing  term,  viz.. 


d2x. 

2. 

dt2 


T(xj.i  - 2xj  * xm>- 


(3.1) 


23  22 

As  N-x*>  Morse  and  Ingard  give  the  Schroedinger  solution  in  our 

notation  as 


23.  P.M.  Morse  and  K.U.  Ingard,  Theoretical  Acoustics  (McGraw-Hill, 
New  York,  1968).  Chap.  3. 
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Figure  2.  Model  for  solving  the  equations  of  motion  for  a one-dimensional,  discrete  lattice. 


x ■ (t)  = y < a J_ . _ (a) 
J L IP  2j-2p'- 

p=_oo  v 


and 
dx . (t)  00 


V»')J2j-2p-dV)  • 

V ° m^p+1  / * 


(3.2) 


at)  “ ( i / 

dt  = l |VpJ2j.2p(wot)t2“oVVl)J2j-2p-l(V)|-  (3’3) 

p = _oo 

where  a^,  are  the  initial  displacement  and  velocity  at  time  t=o, 

J (ui  t)  are  the  Bessel  functions  of  the  first  kind,  of  order  m,  and 
m o / — 

o,  = 2/1  . 

o m 


We  can  guess  from  symmetry  considerations  the  initial  conditions 
which  must  be  imposed  on  the  particles  j < - 1 in  order  that  the  zeroth 
particle  has  the  equation  of  motion  xQ=ut  for  all  time.  Consider  an 

observer  in  a frame  of  reference  moving  at  velocity  u which  is  located 
at  the  origin  of  our  stationary  system  at  time  t=o.  If  the  initial 
conditions  of  the  particles  j > 1 are  the  mirror  image  of  the  particles 
j < - 1 in  this  moving  system,  there  will  be  no  force  on  the  zeroth 
particle.  Therefore,  we  assume  the  following  initial  conditions. 


2u-V 


a =0,  V =u 
o o 


(3.4) 


It  can  be  shown  that 

00  00 
l Vm  = 2uIpI  - u + l Vm  * p < 0 

m=p+l  m= | p | 


(3.5) 


Rearrangement  of  Eq.  (3.2)  using  Eq.  (3.4)  and  Eq.  (3.5)  results  in  the 

solution  to  the  semi-infinite  chain  with  the  boundary  condition 

x =ut,  viz. , 
o 


x.  (t) 


2 |ap(J2j-2p(%t)  - J2j  + 2pC“ot)) 

p=l 

+ — / ? v \ / Jr  , , (w  t)  + J-.  , . (to  t)\ 

w I mil  2j  + 2p-l  o J 2j-2p+l  o J 

°ym=P  / v 

. ^ (2p-l)  J2j,2p.,(»0t)  j ■ > >“ 


(3.6) 
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and 


dx. (t) 
dt 


0)  “ / 

= J (a  -a  .)  (J..  . ,(u  t)  + J_.  « . (u>  t) 

2 L p p+1  l 2j-2p-l  o 2j+2p+l  o 

p=o 

* J (Vp  J2j-2p<“ot>  * <2u-Vp>  J2j*2p(“o«) 
p=o  X 


) 


- u (lOQt)  , j > o 

where  the  following  relationship  has  been  used 

at  J.V  - Vi«V>)  ■ 

Other  useful  relationships  of  Bessel  functions  are 

J-m(wot)  = m inteier 

00  00 

i = y j_  (u)  t)  = j (u  t)  + 2 y j0  (oj  t)  , 

L 2m v o ' ov  o ' L 2mv  o ' 
m=-°°  m=l 

00 

V = 2 l (2m+l)  J2m+1(0)ot)  , 

m=o 

and 


J (o)  = 6 
nr  om 

With  the  above  relationships  one  can  easily  verify  that  x.(o) 

x . (o)  = V.. 

J J 

From  Eq.  (3.6)  and  Eq.  (3.9c)  one  can  show  that 

xo(t)  = zr  E (2p+1)  J2P+i(u,ot)  = ut  ’ 
o p=o  * 

which  is  the  correct  boundary  condition. 


(3.7) 

(3.8) 

(3.9a) 

(3.9b) 

(3.9c) 

(3 . 9d) 

V 

(3.10) 
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A special  solution  of  Eq.  (3.6)  used  in  this  investigation  for  the 
quiescent  lattice,  (a^=o,  y.=o  for  j > o,  Vq=  u) , is 


1 jW  I J2j*2p.l("ot:i  ' 

o p=o  r 


(3.11) 


Another  special  solution  used  in  this  investigation  of  the  semi- 

ut 

infinite  chain  with  the  boundary  condition  x = ut  - -x  , t > t 

o 2 wax 

can  be  determined  by  making  the  following  observations.  Let 

ut 

Yj(t)  = x.(t)  - — p , (3.12) 

where  t = t-t  so  that  the  initial  conditions  at  t=o  are 
max 

Y0(o)  = o,  Yo(o)=u,  Yj(o)=a.  - -^1  , 


and 


(3.13) 


Yj  (o)  = ^j , (j  > o)  , 


+■  Vi 

where  and  ^ are  the  position  and  velocity  of  the  j n particle  at 
t=o. 


Yj (t)  satisfies  an  equation  of  the  same  form  as  Eq.  (3.1),  and  has  the 

same  initial  condition|  at  t=o  as  x.(o)  in  Eq.  (3.6).  Therefore,  the 

max  . ^ 


x.  (t-t  ) = Y 
j max  L.  , 
P=1 


^ 2 x ° 

(/  Ut  \ 

m 

Vp"  2 ) 

J-.  - (u>  (t-t  )' 

2j-2p\  o max  > 

l - J-.  _ (u  (t-t  )') 

2j+2p\  ov  max'/ 

' / r- 

- 

+ — ( I V \ 

o\m=p  J 

+ J_.  ~ . (m  (t-t  )) 

2j -2p+l V o max'/ 


J2j  + 2p-l(wo*‘t_tmax0 


2u 


(3.14) 


ut 


* iT  f2'-1'  J2j.2p-l(“o(‘-W)|  * -F  » < 


max 
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A. 


and 


dx.(t-t  ) to  oo 
j max  " - 


dt 


2 p=o  ('aj"aj  + 1')  ^ J2j-2p-l(a)o(t_tmax^ 

+2p+l  ^“o^t-tmax^  j + jj=o  | ^p  J2j -2p^>o(t"tmax)) 

+ (2u’V  J2j+2p(%(t'W)|-  u J2j  ^o^w)  * 


+ J2j 


t > t 


max 


(3.15) 


The  above  solutions  are  matched  at  t=t  to  the  solutions  which  have 

max 


a different  boundary  condition  between  t=o  and  t=t 


max 

Let  us  now  consider  a special  solution  to  Eq.  (3.1)  with  the  follow- 
ing initial  conditions  at  t=o,  viz.,  aj=0,  vq=0,  v y | are  arbitraty, 

and  v p | are  to  be  determined,  with  the  boundary  condition 


and 


. ut 
ut  max 

o 2 2 it 


dx 

o u u 
dt  - 2 “ 2 C0S 


■(e)- 

fc) 


0 < t < t 


max 


(3.16) 


(3.17) 


Under  these  conditions  Eq.  (3.3)  and  Eq.  (3.17)  reduces  to 


(r-)V  \ 

I \ max  / / p=l 


VV  J2p(V>  ’ 


(3.18) 


where  the  V are  to  be  determined  for  arbitrary  V . The  cos 
-P  P 

24 

has  an  expansion  in  terms  of  Bessel  functions  depending  on  the  range 


(rM 

\ max  / 


t u> 
max  o 


, viz.. 


24.  Handbook  of  Mathematical  Functions,  edited  by  M.  Abramowitz  and 
I.  Stegun  (Nat’l  Bur.  Std.,  WASH,  DC,  1964),  Chap.  9). 
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(3.19) 


cos(u)  t sinB)  = J ( a)  t J + 2 7 J_  (u>  t)  cos  2n  6, 

o oo  L . 2n  o 

n=l 


for  the  case 


t 0) 

max  o 


< 1 where  sin  3 = 


t 0) 
max  o 


and 


cos 


(wnt  cosha)  = Jo(a)Qt)  +2  £ (-1)  J7v(u)nt)  cosh  2k  a,  (3.20) 


k=l 


2k  ^ o 


for  the  case 


t w 
max  o 


> 1 where  cosh  a = 


t u 
max  o 


With  the  aid  of 


Eq.  (3.9b)  and  the  above  equations,  Eq.  (3.18)  can  be  written  as 

OO 

u I 


y (v  + v ) j,  (u  t)  = 

p=i  p -p  2p  ° 


p=1  (1 “Cos  2pB)  J2p(wot) 


u l (l-(-l)p  cosh  2po)  J (a)  t) 


p=l 


2p v o 


(3.21) 


dxQ 

is  an  example  of  an  entire  function.  Therefore,  it  has  a unique 

Neumann's  expansion.  We  can  conclude  that  the  relationship  between  the 

arbitrary  V and  the  determined  V is 
P -p 


u(l-cos2pB) , 


V + V 
P ~P 


t u) 
max  o 


< 1 


(3.22) 


u(l-(-l)H  cosh  2pa) , > 1 

t u 
max  o 


A note  of  caution  is  appropriate  at  this  point.  If  one  tries  to  oLtain 
Eq.  (3.22)  by  multiplying  Eq.  (3.18)  by 


J_  / (u  t) 
2p  v o J 


u>  t 
o 


and  then  integrating  from  0 to  ® with  the  aid  of  Kapteyn's  orthogonality 

. 25 

relation  , 


25.  W.  Kapteyn,  "Sur  Quelques  Integrales  Definies  Contenant  Des  Fonctions 
De  Bessel",  Archives  Neerlandaises  Des  Sciences  Exactes  et  Naturelles, 
VI,  103  (1901). 
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dt 
o 

a different  result  is  obtained.  This  procedure  requires  a questionable 
interchange  of  two  limit  operations,  one  being  the  infinite  series,  the 
other  the  infinite  (improper)  integral.  Apparently  Eq.  (3.20)  is  not 

26 

uniformly  convergent  as  pointed  out  by  Gautschi  . The  following 
summations  are  used  for  the  final  solution,  viz.. 


f 


-xt 


Jnct) 


■ ir(  - x)  • " 


> 1 


(3.23) 


P 

l 

m=l 


cos (2mB) 


sin[  (2p+l)8  ] _ 1_ 

2 sin  6 " 2 


(3.24) 


and 


P 

l (-1)1"  cosh(2ma) 
m=l 


(~l)Pcosh(2p+l)ct  - cosha 
2 cosh  a 


(3.25) 


Therefore,  the  solution  to  Eq.  (3.1)  for  a semi-infinite  chain  with  the 
boundary  condition  Eq.  (3.16),  and  the  initial  conditions 


aj=0,  V | j 1=0  ! 


(3.26) 


is 


l J2jt2pfl(uot]<(pty- 
o p-o  r 


sin[(2p+l)3]  tt  ^ . 

2sin6  ’t  a) 
max  o 


(-l)Fcosh(2p+l)a 

2 cosh  a ’ t 


u 

max  o 


> 1 


(3.27) 


and 


dx . (t) 
J 

dt 


l 

p=l 


(l-cos2pg)  , - - < 1 

max  ^ o 

(l-(-l)Pcosh2pa),  > 1 

max  o 


J (o)  t)  . (3.28) 
2j+2p'  o 


26.  W.  Gautschi,  private  communication. 
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Eq.  (3.27)  and  Eq.  (3.28)  at  t=t  become  the  initial  conditions 
a.,  I/,  for  Eq.  (3.14)  and  Eq.  (3. 15)  ,m§iz. , 


and 


a.  = x. (t  ) , 
j j v max 


dx . (t  ) 
y _ j max 

j 


dt 


(3.29) 


(3.30) 


We  now  have  the  matched  solution  for  a semi- infinite  chain  whose  end 
particle  is  accelerated  in  a prescribed  manner  from  zero  to  constant 
velocity  u in  time  t 

max 


3.2  Propagation  in  the  Initially  Quiescent  Lattice. 

If  we  assume  that  initially  all  particles  except  the  first,  which 
travels  at  velocity  u (xQ=ut) , are  at  rest  in  their  equilibrium  positions, 

we  have  Eq.  (3.11).  The  velocity  of  the  particle  is 

dx. (t) 
dt 

The  infinite  series  in  Eq.  (3.11)  and  Eq.  (3.31)  converge  rapidly,  and 
a computer  program  was  written  to  evaluate  the  sums. 

In  Figure  3,  we  have  plotted  the  nondimensionalized  velocity  and 

displacement  of  the  10t^  particle  as  a function  of  nondimensionalized 
time,  i.e., 


u J0 . (<o  t)  + 2u 
2j  o J 


p=l 


J_.  _ (to  t) 
2j+2pv  o 1 


(3.31) 


Sj(t)  = 


Xj  (t) , 


dSj (t) 

dx 


1_ 

u 


dXj (t) 

dt 


T=tO  t 

o 


(3.32) 


The  peak  velocity  of  the  wave  disturbance  arrives  at  the  1C)  particle 
at  a rate  of  one  half  particle  per  unit  of  nondimensional  time.  This 
result  is  the  maximum  normal  - mode  velocity  as  can  be  obtained  from 
the  dispersion  relation  for  the  harmonic  lattice,  viz.. 


io(k)  = ioo  Isin  | (3.33) 

where  u>(k)  is  the  frequency  of  the  normal  mode  with  wave  vector  k. 

The  group  velocity  is 


du>(k)  “0% 
dk  2 


(3.34) 
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Figure  3.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  10tfl 
particle  are  plotted  as  a function  of  time  for  a harmonic 
lattice  with  a rise  time  t =0.0. 


VELOCITY 


where  the  maximum 

du(k)  i 
dk  1 


is  given  by 


(0  fl. 
o o 


(3.35) 


max 

Per  unit  of  nondimensionalized  time  T=wot,  the  speed  is  aQ/2,  or  one 

half  lattice  spacing  per  unit  tau.  At  a later  time  the  wave  has 

reached  the  60^  particle,  as  can  be  seen  in  Fig.  4.  The  dispersive 
nature  of  the  wave  is  evident. 


The  results  for  all  cases  where  t ^0.0  are  obtained  from  Eq.  (3.27) 

max 

and  Eq.  (3.28)  for  t < tmax*  anc*  from  Eq.  (3.14)  and  Eq.  (3.15)  for 

t > t . When  the  rise  time  t =1.0.  we  notice  little  difference  in 
max  max 

Figure  5 from  the  case  t =0  in  Figure  3.  For  a rise  time  t =6.0  we 

max  6 max 

notice  that  each  amplitude  of  the  successive  peaks  in  the  wave  train  at 

the  lO1"^  particle  in  Figure  6 is  less  than  the  corresponding  amplitude 

at  the  10t  particle  in  Figure  3.  When  this  wave  reaches  the  60th 
particle  in  Figure  7,  each  amplitude  is  not  so  great  as  the  corresponding 

amplitude  at  the  bO^particle  in  Figure  4.  This  trend  continues  as  the 
rise  time  increases.  For  a rise  time  t =12.00  each  amplitude  in  the 

t h max 

wave  train  for  the  10  particle  in  Figure  8 is  very  low.  By  the  time 

the  wave  reaches  the  60  particle  in  Figure  9 its  amplitudes  are  much 

smaller  than  at  the  60^  particle  in  Figure  4.  When  the  wave  reaches 
t h 

the  140  particle  in  Figure  10,  it  still  does  not  have  the  amplitudes 
of  the  OO1*1  particle  in  Figure  4,  although  its  amplitudes  have  increased 
from  their  value  at  the  60^  particle  in  Figure  9. 

Now  we  can  make  some  general  observations  about  the  velocity-time 
trajectories  at  specific  particles  near  the  surface  in  the  harmonic 
approximation.  Each  amplitude  of  the  successive  peaks  in  the  wave  train 
is  less  than  the  corresponding  amplitude  for  the  instantaneous  com- 
pression case  by  an  amount  which  is  inversely  proportional  to  the  rise 
time  This  phenomenon  most  likely  results  from  the  fact  that  the 

total  computed  work  done  on  the  zeroth  particle  at  the  time  when  the 
first  velocity  peak  is  at  a specific  particle  is  less  than  the  corre- 
sponding instantaneous  acceleration  case.  Therefore,  surface  effects 
will  persist  fir  hundreds  of  atoms  if  the  rise  times  are  large  enough. 

It  should  be  pointed  out  that  wave  propagation  in  a quiescent  harmonic 
lattice  has  more  applicability  to  propagation  on  a thermal  background 

than  one  might  at  first  suspect.  Powell  and  Batteh  have  concluded 
that  in  a harmonic  lattice  the  ensemble  average  taken  over  many  initial 
conditions  will  lead  to  the  same  results  for  the  average  velocity  and 
displacement  as  for  the  case  in  which  the  initial  conditions  are  zero. 
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Figure  4.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  60^ 

particle  are  plotted  as  a function  of  time  for  a harmonic 

lattice  with  a rise  time  t =0.0. 

max 
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Figure  6.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  10*^ 

particle  are  plotted  as  a function  of  time  for  a harmonic 

lattice  with  a rise  time  t =6.0. 

max 
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Figure  7.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  60th 

particle  are  plotted  as  a function  of  time  for  a harmonic 

lattice  with  a rise  time  t =6.0. 

max 
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Figure  8.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  10th 
particle  are  plotted  as  a function  of  time  for  a harmonic 
lattice  with  a rise  time  t =12.00. 
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Figure  10.  In  dimensionless  units  the  velocity  which  oscillates 
about  unity  and  the  displacement  from  equilibrium  of 
til 

the  140  particle  are  plotted  as  a function  of  time 

for  a harmonic  lattice  with  a rise  time  t =12.00. 

max 
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4.  NONDIMENSIONALIZED  EQUATIONS 
FOR  ANHARMONIC  CASE  AND  METHOD  FOR  SOLUTION 


If  the  Morse  potential  of  Eq.  (2.9)  is  used  in  Newton's  second 


law  of  Eq.  (2.3),  and  F.  is  set  to  zero,  we  can  obtain  the  nondi- 


mensionalized  equationsJof  motion,  viz., 

‘^[exp{'Vsj+rsj)}“  exp{*  2\rsj+rsj^ } 

- exp{-Am(sj-s._1)}  + exp{-  2AJs.-s._j)} 


d2s.  (T) 

^7" 


(4.1) 


.3  > i 


where  the  definitions  in  Eq.  (3.32)  are  used  as  well  as 

2 


A = au-  , D = 


moi 


m w 


8a 


(4.2) 


The  boundary  condition  for  the  zeroth  particle  can  be  chosen  as  desired. 
In  this  paper  three  cases  are  investigated. 


So  = T»  t > ° 


S = T - 1/2  T , T > T 
o max  max 


0 < t < r 


max 


s = 


o 2t 


, (K  t < t ;s  = t - 1/2  t , t > T 
« 1 max 


max 


max  o 


max 


(4.3a) 


(4.3b) 


(4.3c) 


where  t = u t 

max  o max 


Eq-  (4j-Land  E<?‘  constitute  a set  of  N coupled,  nonlinear,  second- 

order  differential  equations  which  must  be  solved  numerically  for  various 
values  of  the  parameters  under  consideration.  These  equations  can  be 
converted  to  a set  of  2n  first  order  equations,  viz.. 


J J 


*i  ‘ «m  | 'XP  [-2A»<sj-sj-l)]-  exp[-Am(sj-5j-l>] 

[-2Am<sj.rsj)J*  “p[-Vsj.i-sj>]}’j » 


- exp 
d2s 


(4.4) 


Vo  = 


dx 
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where  the  dot  denotes  differentiation  with  respect  to  the  dimensionless 
t ime  t . 

To  solve  Eq.  (4.4)  we  employed  a computer  program  developed  by 

Powell  and  Batteh13,  modified  to  take  into  account  the  different 
boundary  conditions  in  Eq.  (4.3).  This  program  uses  a fourth-order 
27 

Runge-Kutta  scheme  . Given  the  values  of  the  functions  on  the  left- 
hand  side  of  Eq.  (4.4)  at  time  T,  this  method  approximates  their  values 
at  time  x+  At  by  a fourth-order  polynominal  in  At. 

The  harmonic  limit  of  Eq.  (4.4)  occurs  as  Am  tends  to  zero.  When 

= 0.0001,  good  agreement  with  the  harmonic  case  was  obtained.  The 

interested  reader  can  obtain  more  details  on  the  program  from  the  above 
reference  13. 

The  remainder  of  this  paper  will  discuss  the  results  of  the  numeri- 
cal solution  of  Eq.  (4.4)  for  different  cases. 


5.  PROPAGATION  IN  THE  INITIALLY 
QUIESCENT,  ANHARMONIC  LATTICE 

In  this  section  we  will  discuss  the  numerical  solution  of  Eq.  (4.4) 

for  various  values  of  the  anharmonicity  parameter  A and  rise  time  t 

J r m max 

For  this  discussion  the  following  definitions  are  needed.  The  term 
solitary  wave  means  a localized  traveling  wave  of  constant  shape  and 
amplitude.  In  this  report,  the  term  soliton  describes  a nonlinear 
solitary  wave  which  emerges  from  a collision  with  a similar  pulse, 
retaining  the  same  shape  and  speed  it  had  initially.  Finally,  the  term 
envelope  soliton  describes  an  envelope  of  constant  speed  imposed  on  a 
solitary  wave  train  with  its  own  carrier  speed. 

The  results  are  divided  into  four  parts,  each  of  which  represents  a 
different  boundary  condition  on  the  zeroth  particle.  When  the  zeroth 
particle  travels  at  constant  velocity,  we  observe  a solitary  wave  train 
at  the  head  of  the  velocity  trajectories  of  individual  particles  and  an 
oscillatory  tail  which  persists  in  the  long-time  limit.  A sinusoidal  accel- 
eration on  the  zeroth  particle  produces  what  may  be  an  envelope  soliton 
traveling  much  slower  than  the  shock  wave.  The  same  behavior  is  observed 
for  a ramp  acceleration.  Finally,  when  the  zeroth  particle  is  decelerated 

to  zero,  an  example  of  solitons  spreading  out  in  time  is  observed. 


27.  B.  Carnahan,  H.A.  Luther,  and  J.O.  Wilkes,  Applied  Numerical 
Methods  (W  y,  New  York,  1969),  Chap.  6. 
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5 . 1 Zeroth  Particle  Travels  at  Constant  Velocity 


Let  us  look  at  the  results  of  solving  Eq.  (4.4)  under  the  boundary 
condition  Eq.  (4.3a).  This  physically  corresponds  to  the  case  where  the 
zeroth  particle  is  pushed  at  constant  speed  for  all  time.  For  the 
anharmonicity  parameter  Am=0.2,  Tmax=0>  an^  At=0.05  the  velocity  and 

position  of  selected  particles  is  plotted  as  a function  of  tau.  As  the 

wave  travels  from  the  10th  particle  to  the  90th  particle  in  Figure  11 

to  Figure  14,  we  see  the  amplitudes  of  the  successive  velocity  peaks 

near  the  head  of  the  shock  wave  gradually  develop  into  a solitary  wave 

9 13  13 

train  as  was  pointed  out  by  Tasi  and  Powell  and  Batteh  . The  latter 

have  plotted  the  maximum  particle  velocity  behind  the  shock  front  for 
A =0.2  and  A =1.0  at  a time  when  the  front  is  approximately  located  at 
m th  m 

the  480  particle.  Their  result  indicates  that  the  amplitudes  of  the 
leading  solitary  waves  approach  a dimensionless  velocity  of  2.0.  The 
leading  amplitudes  for  the  case  A =1.0  approaze  the  value  of  2.0  much 

sooner  than  the  case  Am=0.2.  The  oscillatory  tail  is  evident  in  Figure 

11  and  Figure  12  as  was  noticed  by  Tasi^.  A similar  oscillatory  tail 
was  reported  by  Zabusky  in  his  numerical  solution  to  the  Korteweg-de 
28 

Vries  equation 

For  the  anharmonicity  parameter  Am=1.0,  Tmax=0. 0, At=0. 025,  and 

boundary  condition  Eq.  (4.3a),  we  see  that  the  solitary  wave  train  at 

the  head  of  the  shock  wave  has  nearly  formed  at  the  20th  particle 

in  Figure  15.  As  the  wave  passes  the  80th  and  240th  particle  in 
Figure  16  and  Figure  17  we  notice  that  the  leading  amplitudes  of  sue- 

th 

cessive  velocity  peaks  have  increased  slightly  over  the  20  particle. 

For  all  these  cases  we  notice  that  the  amplitudes  of  the  oscillatory 
tail  at  times  long  after  the  shock  has  passed  has  not  approached  zero 

as  was  the  case  for  the  10th  and  20th  particle  for  A =0.2.  A careful 

th  th  m 

examination  of  the  80  n and  240L  velocity  trajectories  will  also  show 
a slight  spreading  of  the  leading  solitary  waves.  The  greater  the 
amplitude  the  faster  they  travel.  This  property  is  a characteristic  of 
solitons.  It  was  decided  to  observe  the  oscillatory  tail  at  times  long 

after  the  shock  wave  has  passed.  The  lSt  and  5th  particles  were  observed 

for  times  greater  than  200  in  Figure  18  and  Figure  19.  The  maximum 

amplitude  of  the  1st  particle  is  approximately  1.2  while  that  of  the  5th 

particle  is  approximately  1.5.  The  oscillatory  tail  for  the  20  ' particle 


28.  N.J.  Zabusky,  "Solitons  and  Bound  States  of  the  Time-Independent 
Schrodinger  Equation",  Phys.  Rev.  168 , 124  (1968). 
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Figure  11.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  10t*1 
particle  are  plotted  as  a function  of  time  for  a lattice 
with  an  anharmonicity  parameter  Am=.2  and  a rise  time 
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Figure  12.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  20t*1 
particle  are  plotted  as  a function  of  time  for  a lattice 
with  an  anharmonicity  parameter  A =.2  and  a rise  time 
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Figure  13.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  60*^ 
particle  are  plotted  as  a function  of  time  for  a lattice 
with  an  anharmonicity  parameter  A =.2  and  a rise  time 

A 
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Figure  14.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  90*  ^ 

particle  are  plotted  as  a function  of  time  for  a lattice 

with  an  anharmonicity  parameter  A =.2  and  a rise  time 

t =0.0.  m 

max 
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Figure  15.  In  dimensionless  units  the  velocity  which  oscillates  about 

t h 

unity  and  the  displacement  from  equilibrium  of  the  20 
particle  are  plotted  as  a function  of  time  for  a lattice 
with  an  anharmonicity  parameter  Am=1.0  and  a rise  time 
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Figure  18.  In  dimensionless  units  the  velocity  which  oscillates  about 

cf 

unity  and  the  displacement  from  equilibrium  of  the  1 
particle  are  plotted  as  a function  of  time  for  a lattice 
with  an  anharmonicity  parameter  A =1.0  and  a rise  time 
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Figure  19.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  5^ 
particle  are  plotted  as  a function  of  time  for  a lattice 


A =1.0  and  a rise  time 
m 
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in  Figure  20  was  observed  for  a much  longer  period  of  time,  and  its 

til 

maximum  amplitude  was  found  to  be  slightly  greater  than  the  5 
particle. 

5.2  Zeroth  Particle  Accelerated  Sinusoidally 

Let  us  compare  the  results  of  solving  Eq.  (4.4)  under  the  boundary 
condition  Eq.  (4.3b)  to  the  ones  obtained  under  the  boundary  condition 
Eq.  (4.3a).  This  physically  corresponds  to  the  case  where  the  zeroth 
particle  is  accelerated  sinusoidally  from  zero  velocity  to  its  final 
velocity  in  time  t . For  the  anharmonicity  parameter  Am=0.2, 

x =12.0,  and  At=0.5,  we  observe  that  the  amplitude  of  the  velocity 
msx 

of  the  shock  wave  at  the  20LI  particle  in  Figure  21  is  less  than  the 

t h 

case  t =0.0  in  Figure  12.  As  the  shock  wave  reaches  the  60 L 1 particle 
max  6 v 

the  leading  peak  in  Figure  22  is  the  same  as  Figure  13,  but  the  succeed- 
ing pulses  are  less  in  ampl^ude.  According  to  our  results,  by  the  time 

the  wave  has  reached  the  80  n particle,  the  cases  t =0.0  and  x =12.0 

r max  max 

look  the  same.  Therefore,  for  a low  anharmonicity  parameter  such  as 

Am=0.2  each  amplitude  of  the  successive  peaks  in  the  wave  train  for 

surface  particles  is  less  than  the  corresponding  amplitude  for  the 
instantaneous  compression  case  by  an  amount  which  is  inversely  propor- 
tional to  the  rise  time  rm  . This  phenomenon  is  analogous  to  the 

harmonic  case.  The  most  distinguishing  factor  between  the  harmonic 
and  slightly  anharmonic  cases  is  the  steeper  and  narrower  pulse  width 
for  the  latter,  which  gradually  develops  into  a solitary  wave  train. 

For  the  anharmonicity  parameter  Am=1.0,  Tmax=12.0  and  At=.025  not 

only  do  we  see  the  characteristics  already  described  at  lower  anharmoni- 
cities,  but  also  the  appearance  of  a dip  in  the  amplitude  of  the 
velocity  at  a time  after  the  shock  wave  has  passed  an  individual  parti- 
cle. Perhaps,  the  disturbance  is  an  envelope  soliton.  In  Figure  23 

to  Figure  27  we  observe  the  envelope  traveling  past  the  20^,  30^,  40t  , 

and  60*^  particle  at  approximately  one  sixth  the  speed  of  the  shock 

front.  The  amount  of  computer  time  prevented  us  from  fol^wing  the 

envelope  past  325  tau.  The  velocity  trajectory  of  the  60  particle 

has  the  same  profile  for  r =0  as  for  t =12.0  except  in  the  vicinity 
r max  max  r 

of  the  envelope  from  the  time  the  shock  wave  arrives  until  325  tau. 

The  profiles  of  particles  greater  than  60  look  alike  up  to  325  tau  as 

the  envelope  would  be  expected  to  appear  at  a later  time.^In  Figure 

28  we  show  the  initial  solitary  waves  arriving  at  the  380  particle 

for  the  case  t =12.0,  but  we  expect  the  same  initial  profile  to  hold 
max  r r 

for  other  values  of  tma  . The  380t  particle  is  the  farthest  point  in 
the  lattice  for  which  we  were  able  to  compute  a meaningful  profile. 
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Figure  21.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  20t*1 
particle  are  plotted  as  a function  of  time  for  a lattice 
with  an  anharmonicity  parameter  A -.2  and  a rise  time 
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Figure  22.  In  dimensionless  units  the  velocity  which  oscillates  about 

x.  L 

unity  and  the  displacement  from  equilibrium  of  the  6Cr 
particle  are  plotted  as  a function  of  time  for  a lattice 
with  an  anharmonicity  parameter  A =.2  and  a rise  time 

- *1  n ^ 
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Figure  23. 


In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  20^ 

particle  are  plotted  as  a function  of  time  for  a lattice 

x/ith  an  anharmonicity  parameter  A =1.0  and  a rise  time 

r =12.0.  m 

max 
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We  wanted  to  verify  that  the  envelope  would  occur  for  other 
anharmonicities,  and  at  the  same  time  investigate  the  effect  of  different 
rise  times  xmax.  For  the  anharmonicity  parameter  Am=1.2,  and  At=0.025 

we  looked  at  the  ninth  particle  for  different  values  of  Tmax-  For 

t =4.0  we  see  no  envelope  at  all  in  Figure  29.  This  result  is  most 

likely  explained  by  the  fact  that  as  ^max  approaches  zero  no  envelope 

has  been  reported  by  any  investigator.  For  values  of  Tmax=8.0,  12.0,  and 

20.0  in  Figure  30  through  Figure  32,  respectively,  we  observe  an  envelope 
which  appears  to  deepen  and  appear  at  a later  time  relative  to  the 
arrival  of  the  initial  disturbance  at  the  ninth  particle. 

5 . 3 Zeroth  Particle  Given  a Ramp  Acceleration 

We  wanted  to  make  sure  that  what  seemed  to  be  an  envelope  soliton 
was  not  caused  by  our  choice  of  the  sinusoidal  acceleration  given  to 
the  zeroth  particle,  rather  than  some  other  acceleration,  like  the  ramp 
acceleration  in  Eq.  (4.3c).  For  the  anharmonicity  parameter  Am=1.2, 

and  Ax=.025  an  envelope  was  observed  at  the  ninth  particle  for  Tmax=4-0 

and  Tmax=12.0  in  Figure  33  and  Figure  34,  respectively.  Even  for 

Tmax=^’°  an  enveloPe  was  formed  in  contrast  to  the  case  in  Figure  29 

for  the  same  rise  time.  For  x =12.0  the  envelope  was  deeper  than  the 

max  r r 

case  x =4.0  and  occurred  at  a later  time  after  the  arrival  of  the 
max 

shock  wave. 

5.4  Zeroth  Particle  Decelerated  to  Zero 


Finally,  we  wanted  to  see  the  effect  of  starting  the  zeroth  parti- 
cle at  constant  velocity  and  then  decelerating  it  to  zero  in  time 
Tmax'  ^he  boundary  condition  is  Eq.  (4.3a)  minus  Eq.  (4.3b).  For  the 

anharmonicity  parameter  Am=1.0,  Ax=.025,  and  Tmax=12.0  we  have  a beauti- 
ful example  of  solitons  spreading  out  in  time  starting  at  the  tenth 
particle  in  Figure  35  and  going  to  the  fortieth  particle  in  Figure  36. 


6.  DISCUSSION 

In  this  report  we  have  investigated  the  effect  of  accelerating  the 
end  particle  of  a one-dimensional  lattice  to  its  final  velocity  in  a 
characteristic  time  Tmax-  The  purpose  was  to  determine  the  shock  pro- 
file caused  in  this  manner  and  compare  it  with  the  instantaneous  com- 
pression case.  For  harmonic  lattices  we  noticed  that  in  the  surface 
atoms  each  amplitude  of  the  successive  peaks  in  the  wave  train  was  less 
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Figure  29.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  9^ 
particle  are  plotted  as  a function  of  time  for  a lattice 
with  an  anharmonicity  parameter  A =1.2  and  rise  time 
4. 
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Figure  30.  In  dimensionless  units  the  velocity  which  oscillates  about 

t h 

unity  and  the  displacement  from  equilibrium  of  the  9l 

particle  are  plotted  as  a function  of  time  for  a lattice 

with  an  anharmonicity  parameter  A =1.2  and  rise  time 

t =8.0.  m 

max 
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Figure  31.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  9tfl 

particle  are  plotted  as  a function  of  time  for  a lattice 

with  an  anharmonicity  parameter  A =1.2  and  rise  time 

t =12.0.  m 

max 
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Figure  32.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  9t^1 
particle  are  plotted  as  a function  of  time  for  a lattice 
with  an  anharmonicity  parameter  A =1.2  and  rise  time 
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Figure  33.  In  dimensionless  units  the  velocity  which  oscillates  about 

unity  and  the  displacement  from  equilibrium  of  the  9t*1 
particle  are  plotted  as  a function  of  time  for  a lattice 
with  an  anharmonicity  parameter  A =1.2  and  rise  time 
T =4.0.  A ramp  acceleration  of™ zeroth  particle  is 


VELOCITY 


DISTANCE 


0 40  80  120  160  200 

r 

Figure  34.  In  dimensionless  units  the  velocity  which  oscillates  about 

t h 

unity  and  the  displacement  from  equilibrium  of  the  9L 
particle  are  plotted  as  a function  of  time  for  lattice 
with  an  anharmonicity  parameter  Am=1.2  and  rise  time 

TmaX”12-0'  A ramp  acceleration  of  zeroth  particle  is  used. 
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Figure  35.  In  dimensionless  units  the  velocity  which  propagates  as 

four  solitary  waves  and  the  displacement  from  equilibrium 

of  the  10*^  particle  are  plotted  as  a function  of  time  for 
a lattice  with  an  anharmonicity  parameter  Am=1.0  and 

deceleration  time  Tmax=12.0.  The  zeroth  particle  starts 

out  with  initial  velocity  of  unity  and  is  decelerated  to 
zero. 
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than  the  corresponding  amplitude  for  the  instantaneous  compression  case 
by  an  amount  which  is  inversely  proportional  to  the  rise  time  Tmax-  As 

the  lattice  becomes  slightly  anharmonic,  each  amplitude  of  the  succes- 
sive peaks  in  the  wave  train  passing  through  the  surface  particles 
behaves  in  a manner  which  is  analogous  to  the  harmonic  case.  The  most 
distinguishing  factor  between  the  harmonic  and  slightly  anharmonic 
cases  is  the  steeper  and  narrower  pulse  width  for  the  latter,  which 
gradually  develops  into  a solitary  wave  train  at  the  head  of  the  shock 
wave  while  an  oscillatory  tail  persists.  Finally,  for  an  anharmonicity 
parameter  one  or  greater  we  noticed,  in  addition  to  the  above  anharmonic 
effects,  what  appears  to  be  an  envelope  soliton  propagating  at  a slower 
speed  than  the  shock  wave  for  certain  rise  times. 

As  was  mentioned  in  the  introduction,  a very  important  and  current 
research  area  is  the  process  by  which  the  lattice  returns  to  thermal 
equilibrium,  if  at  all,  after  it  has  been  perturbed  by  a shock  wave 
propagating  on  a thermal  background.  All  the  investigators  of  this 
problem  sample  their  data  behind  the  shock  front  to  determine  whether 
or  not  the  criteria  for  thermal  equilibrium  have  been  satisfied.  The 
question  of  whether  chemical  reactions  in  reactive  materials  occur  in 
the  equilibrated  or  the  non-equilibrated  region  of  the  crystal  is  very 
important  for  detonation  theory.  Surface  chemistry  will  also  be 
influenced  by  the  strength  of  the  disturbance  passing  through  this 
region.  It  is  possible  that  different  rise  times  could  affect  the 

above  phenomena.  Tsai16  has  already  reported  that  instantaneous 
compression  of  the  end-most  particle  causes  a large  increase  in  the 
kinetic  energy  of  the  surface  atoms  which  takes  a long  time  to  ther- 
malize. 

There  also  remains  the  problem  of  why  the  principal  investi- 
gators interpret  their  results  differently.  For  instance,  Tsai14  16 
uses  non-zero  rise  times  and  claims  that  the  shock  profile  is  not 
steady  in  time  for  a three-dimensional  lattice.  He  says  that  there 
is  an  ever-expanding  region  of  non-equilibrium  between  the  shock  front 
and  the  thermally  equilibrated  region  behind  the  front.  On  the  other 
1-2 

hand  Paskin  uses  a zero  rise  time  and  claims  that  the  shock  profile 

is  steady  in  time.  Powell  and  Batteh1^  use  a zero  rise  time  and  have 
found  that  the  shock  profile  is  not  steady  in  a one-dimensional  lattice. 
Unfortunately,  we  did  not  propagate  our  shock  wave  on  a thermal  back- 
ground in  this  investigation.  However,  we  feel  that  the  different 
wave  characteristics  in  the  surface  atoms,  and  the  possibility  of 
envelope  solitons  propagating  at  a slower  speed  than  the  shock  front, 
would  exist,  but  would  be  masked  by  the  thermal  background.  Therefore, 
future  investigators  should  be  aware  of  these  results  as  they  strive  for 
a uniform  interpretation  of  the  shock  wave  problem. 
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